Análisis de expresión diferencial

 

# Load packages 
library(pacman)
pacman::p_load(glue)
pacman::p_load(DT)
pacman::p_load(ggpubr)
pacman::p_load(flextable)
pacman::p_load(officer)
pacman::p_load(dplyr)
pacman::p_load(magick)

# Functions
source("C:/home/rmoldovan/T2D-Meta-Analysis/Functions/Functions.R")
Dir = "/home/rmoldovan/T2D-Meta-Analysis/Data/04DE" 

 

# Plot output
RX_grid_img <- function(paths){
  grid = NULL
  for (i in seq(from = 1, to=length(paths), by=2)){
    a = image = image_read(paths[i])
    b = image_read(paths[i+1])
    row = image_append(c(a, b))
    print(i)
    if(is.null(grid)){
      grid = row
    }else{
      grid = image_append(c(grid,row),  stack = TRUE)
    } }
    return(grid)
    }
# Output table
RX_table <- function(Data,
                     Color1 = "#007FA6",
                     Color2 = "#007FA6" ,
                     Color.l = "black" ,
                     Color.l2 = "#878787",
                     Color.up = "#CB326D",
                     Color.down = "#00AFBB",
                     footer = NULL,
                     font = "Calibri"){
    
  cols=colnames(Data)[-c(1,2)] # Columns with info
  col_keys = colnames(Data) # Nombres unicos de las columnas
  # Header
  header = data.frame(col_keys = col_keys,
                      H1 = c("Contraste", "Contraste", rep("Estudios", length(cols))),
                      H2 = c("Contraste", "Contraste", cols)
                      )
  
  ## Flextable
  # Set info
  ft <- flextable(Data, col_keys = colnames(Data))
  # Set header
  ft <- set_header_df(ft, mapping = header, key = "col_keys" )
  # Format number
  ft <- colformat_num(ft, big.mark = ".", decimal.mark = "," )
  # Set footer
  if(! is.null(footer)){
    ft <- add_footer(ft, values = footer) %>% 
      merge_h(part = "footer")  %>% 
      bold(j=1, part = "footer") %>%
      hline(border = fp_border(width = 2, color = Color1), part = "footer")
  }
  
  # Merge headers
  ft <- merge_h(ft, part = "header") %>% 
    merge_v( part = "header")
  
  # Merge contrastes
  ft <- merge_v(ft, j = c(1,2), part = "body")
  
  # Horizontal lines
  ft <- hline(ft, i=seq(from=3, to=nrow(Data), by=3),
              border = fp_border(width = 1, color = Color.l2), part="body")
  ft <- hline(ft,border = fp_border(width = 2, color = Color1), part = "header")
  
  # Outer bottom
  ft <- hline_bottom(ft, border = fp_border(width = 2, color = Color1))
  
  # Alineamiento
  ft <- flextable::align(ft, align = "center", part = "header")
  ft <- flextable::align(ft, align = "right", part = "footer")
  ft <- flextable::align(ft, j=1, align = "left", part = "footer")
  
  
  # Header design
  ft <- fontsize(ft, i = 1:2, size = 12, part = "header") %>%
    color(i=1, color = Color2, part = "header") %>%
    color(i=2, color = Color.l2, part = "header") %>%
    bold(i=1, bold = TRUE, part = "header")
  # Rotate names
  #ft <- rotate(ft, i = 2, rotation = "btlr", part = "header", align = "bottom")
  
  # Font
  #font(ft, fontname = font, part = "all")
  
  # Matriz de colores
  colormatrix <- ifelse(Data[, cols] == 0, Color.l2, Color.l )
  # Color
  ft <- color(ft, j=cols, color = colormatrix, part = "body") 
  # Color
  x = Data[, 2]
  colors = case_when(
    x == "Up" ~ Color.up,
    x == "Down" ~ Color.down,
    TRUE ~ Color.l
  )
  ft <- color(ft, j=2, color = colors, part = "body") %>%
    bold(j=2, bold = TRUE, part = "body")
  return(ft)
    
}

 

Group

 

# Read data
Res = get(load(glue("{Dir}/DifferentialExpressionGroup.RData")))
St = names(Res)
Toptabs = sapply(names(Res), 
                 function(est) sapply(names(Res[[est]]), 
                                      function(tis) sapply(Res[[est]][[tis]], 
                                                           function(x) x$TopTab, 
                                                           simplify = FALSE),
                                      simplify = FALSE),
                 simplify = FALSE) 
Plots = sapply(names(Res), 
                 function(est) sapply(names(Res[[est]]), 
                                      function(tis) sapply(Res[[est]][[tis]], 
                                                           function(x) x$Plot, 
                                                           simplify = FALSE),
                                      simplify = FALSE),
                 simplify = FALSE) 

 

SAT

 

Tabla
Toptabs_tis = sapply(Toptabs, 
           function(x) x[names(x) == "SAT"],
           simplify = FALSE, USE.NAMES = TRUE) 

Simp = sapply(Toptabs_tis, function(x) RX_ifelse(is.null(x$SAT[1]), NULL, x$SAT),
              simplify = FALSE) 
Toptabs_tis = Filter(Negate(is.null), Simp)

Toptabs_tis_sig = sapply(Toptabs_tis, 
           function(x) sapply(x, 
                              function(x1) x1[x1$"adj.P.Val" <= 0.05,],
                              simplify = FALSE, USE.NAMES = TRUE),
           simplify = FALSE, USE.NAMES = TRUE) 

Resume_tis = sapply(Toptabs_tis_sig, 
                    function(x) sapply(x, 
                              function(x1){ up = length(which(x1$logFC > 0))
                                            down = length(which(x1$logFC < 0))
                                            return(c(up,down, sum(up,down)))},
                              simplify = FALSE, USE.NAMES = TRUE),
           simplify = FALSE, USE.NAMES = TRUE)

rows = unlist(sapply(Resume_tis, function(x) names(x), USE.NAMES = TRUE, simplify = FALSE))
rows = unique(rows)
# Sep
cols = names(Resume_tis)

DT_tis0 = sapply(cols,
                 function(col) sapply(rows,
                                      function(row) RX_ifelse(is.null(Resume_tis[[col]][[row]]), rep("-", 3), Resume_tis[[col]][[row]]),
                                      USE.NAMES = T),
                 USE.NAMES = T)
DT_tis = data.frame("Contrast" = rep(rows, each = 3),
                    "Direction" = rep(c("Up", "Down", "Total"), times = length(rows)),
                    DT_tis0)
#datatable(DT_tis0, caption = "Resultados expresión diferencial SAT")
foot = sapply(names(Toptabs_tis), function(x) nrow(Toptabs_tis[[x]][[1]]))
val = c("Nº total de genes", "Nº total de genes",foot)
names(val)[1:2] = c("Contrast", "Direction")
RX_table(DT_tis, footer = val) # OUT

Contraste

Estudios

E_MEXP_1425

GSE2508

GSE20950

GSE29718

GSE64567

GSE92405

GSE141432

GSE205668

Ob_IS - Np_IS

Up

719

531

-

0

106

987

617

1632

Down

258

413

-

0

37

1404

304

307

Total

977

944

-

0

143

2391

921

1939

Ob_IS.M - Np_IS.M

Up

616

172

-

0

0

674

2

59

Down

198

124

-

0

0

1047

4

4

Total

814

296

-

0

0

1721

6

63

Ob_IS.F - Np_IS.F

Up

57

338

-

0

18

1097

201

0

Down

57

191

-

0

13

1275

167

0

Total

114

529

-

0

31

2372

368

0

(Ob_IS.M - Np_IS.M) - (Ob_IS.F - Np_IS.F)

Up

0

0

-

0

0

763

0

0

Down

0

1

-

0

0

635

2

0

Total

0

1

-

0

0

1398

2

0

Ob_IR - Ob_IS

Up

-

-

1998

-

-

-

3

-

Down

-

-

4670

-

-

-

0

-

Total

-

-

6668

-

-

-

3

-

Ob_IR.M - Ob_IS.M

Up

-

-

248

-

-

-

0

-

Down

-

-

509

-

-

-

0

-

Total

-

-

757

-

-

-

0

-

Ob_IR.F - Ob_IS.F

Up

-

-

1416

-

-

-

2

-

Down

-

-

2736

-

-

-

0

-

Total

-

-

4152

-

-

-

2

-

(Ob_IR.M - Ob_IS.M) - (Ob_IR.F - Ob_IS.F)

Up

-

-

0

-

-

-

0

-

Down

-

-

0

-

-

-

0

-

Total

-

-

0

-

-

-

0

-

Ob_IR - Np_IS

Up

-

-

-

-

-

-

2171

-

Down

-

-

-

-

-

-

1175

-

Total

-

-

-

-

-

-

3346

-

Ob_IR.M - Np_IS.M

Up

-

-

-

-

-

-

15

-

Down

-

-

-

-

-

-

32

-

Total

-

-

-

-

-

-

47

-

Ob_IR.F - Np_IS.F

Up

-

-

-

-

-

-

1491

-

Down

-

-

-

-

-

-

735

-

Total

-

-

-

-

-

-

2226

-

(Ob_IR.M - Np_IS.M) - (Ob_IR.F - Np_IS.F)

Up

-

-

-

-

-

-

0

-

Down

-

-

-

-

-

-

3

-

Total

-

-

-

-

-

-

3

-

Nº total de genes

21367

19891

21367

19975

16589

21367

35795

35949

 


Plots
E_MEXP_1425
contrasts = c('Ob_IS-Np_IS', 'Ob_IS.M-Np_IS.M', 'Ob_IS.F-Np_IS.F', '(Ob_IS.M-Np_IS.M)-(Ob_IS.F-Np_IS.F)')
paths = glue("{Dir}/Plots/E_MEXP_1425_SAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3

 


GSE2508
contrasts = c('Ob_IS-Np_IS', 'Ob_IS.M-Np_IS.M', 'Ob_IS.F-Np_IS.F', '(Ob_IS.M-Np_IS.M)-(Ob_IS.F-Np_IS.F)')
paths = glue("{Dir}/Plots/GSE2508_SAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3

 


GSE20950
contrasts = c('Ob_IR-Ob_IS', 'Ob_IR.M-Ob_IS.M', 'Ob_IR.F-Ob_IS.F', '(Ob_IR.M-Ob_IS.M)-(Ob_IR.F-Ob_IS.F)')
paths = glue("{Dir}/Plots/GSE20950_SAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3

 


GSE29718
contrasts = c('Ob_IS-Np_IS', 'Ob_IS.M-Np_IS.M', 'Ob_IS.F-Np_IS.F', '(Ob_IS.M-Np_IS.M)-(Ob_IS.F-Np_IS.F)')
paths = glue("{Dir}/Plots/GSE29718_SAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3

 


GSE64567
contrasts = c('Ob_IS-Np_IS', 'Ob_IS.M-Np_IS.M', 'Ob_IS.F-Np_IS.F', '(Ob_IS.M-Np_IS.M)-(Ob_IS.F-Np_IS.F)')
paths = glue("{Dir}/Plots/GSE64567_SAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3

 


GSE92405
contrasts = c('Ob_IS-Np_IS', 'Ob_IS.M-Np_IS.M', 'Ob_IS.F-Np_IS.F', '(Ob_IS.M-Np_IS.M)-(Ob_IS.F-Np_IS.F)')
paths = glue("{Dir}/Plots/GSE92405_SAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3

 


GSE141432
contrasts = c('Ob_IR-Ob_IS', 'Ob_IR-Np_IS', 'Ob_IS-Np_IS', 'Ob_IR.M-Ob_IS.M', 'Ob_IR.M-Np_IS.M', 'Ob_IS.M-Np_IS.M', 'Ob_IR.F-Ob_IS.F', 'Ob_IR.F-Np_IS.F', 'Ob_IS.F-Np_IS.F', '(Ob_IR.M-Ob_IS.M)-(Ob_IR.F-Ob_IS.F)', '(Ob_IR.M-Np_IS.M)-(Ob_IR.F-Np_IS.F)', '(Ob_IS.M-Np_IS.M)-(Ob_IS.F-Np_IS.F)')
paths = glue("{Dir}/Plots/GSE141432_SAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3
## [1] 5
## [1] 7
## [1] 9
## [1] 11

 


GSE205668
contrasts = c('Ob_IS-Np_IS', 'Ob_IS.M-Np_IS.M', 'Ob_IS.F-Np_IS.F', '(Ob_IS.M-Np_IS.M)-(Ob_IS.F-Np_IS.F)')
paths = glue("{Dir}/Plots/GSE205668_SAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3

 


VAT

 

Tabla
Toptabs_tis = sapply(Toptabs, 
           function(x) x[names(x) == "VAT"],
           simplify = FALSE, USE.NAMES = TRUE) 

Simp = sapply(Toptabs_tis, function(x) RX_ifelse(is.null(x$VAT[1]), NULL, x$VAT),
              simplify = FALSE) 
Toptabs_tis = Filter(Negate(is.null), Simp)

Toptabs_tis_sig = sapply(Toptabs_tis, 
           function(x) sapply(x, 
                              function(x1) x1[x1$"adj.P.Val" <= 0.05,],
                              simplify = FALSE, USE.NAMES = TRUE),
           simplify = FALSE, USE.NAMES = TRUE) 

Resume_tis = sapply(Toptabs_tis_sig, 
                    function(x) sapply(x, 
                              function(x1){ up = length(which(x1$logFC > 0))
                                            down = length(which(x1$logFC < 0))
                                            return(c(up,down, sum(up,down)))},
                              simplify = FALSE, USE.NAMES = TRUE),
           simplify = FALSE, USE.NAMES = TRUE)

rows = unlist(sapply(Resume_tis, function(x) names(x), USE.NAMES = TRUE, simplify = FALSE))
rows = unique(rows)
# Sep
cols = names(Resume_tis)

DT_tis0 = sapply(cols,
                 function(col) sapply(rows,
                                      function(row) RX_ifelse(is.null(Resume_tis[[col]][[row]]), rep("-", 3), Resume_tis[[col]][[row]]),
                                      USE.NAMES = T),
                 USE.NAMES = T)
DT_tis = data.frame("Contrast" = rep(rows, each = 3),
                    "Direction" = rep(c("Up", "Down", "Total"), times = length(rows)),
                    DT_tis0)
#datatable(DT_tis0, caption = "Resultados expresión diferencial VAT")
foot = sapply(names(Toptabs_tis), function(x) nrow(Toptabs_tis[[x]][[1]]))
val = c("Nº total de genes", "Nº total de genes",foot)
names(val)[1:2] = c("Contrast", "Direction")
RX_table(DT_tis, footer = val) # OUT

Contraste

Estudios

GSE20950

GSE29718

Ob_IR - Ob_IS

Up

618

-

Down

12687

-

Total

13305

-

Ob_IR.M - Ob_IS.M

Up

229

-

Down

5556

-

Total

5785

-

Ob_IR.F - Ob_IS.F

Up

455

-

Down

9946

-

Total

10401

-

(Ob_IR.M - Ob_IS.M) - (Ob_IR.F - Ob_IS.F)

Up

0

-

Down

0

-

Total

0

-

Ob_IS - Np_IS

Up

-

0

Down

-

2

Total

-

2

Ob_IS.M - Np_IS.M

Up

-

0

Down

-

0

Total

-

0

Ob_IS.F - Np_IS.F

Up

-

0

Down

-

0

Total

-

0

(Ob_IS.M - Np_IS.M) - (Ob_IS.F - Np_IS.F)

Up

-

0

Down

-

0

Total

-

0

Nº total de genes

21367

19975

 


Plots
GSE20950
contrasts = c('Ob_IR-Ob_IS', 'Ob_IR.M-Ob_IS.M', 'Ob_IR.F-Ob_IS.F', '(Ob_IR.M-Ob_IS.M)-(Ob_IR.F-Ob_IS.F)')
paths = glue("{Dir}/Plots/GSE20950_VAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3

 


GSE29718
contrasts = c('Ob_IS-Np_IS', 'Ob_IS.M-Np_IS.M', 'Ob_IS.F-Np_IS.F', '(Ob_IS.M-Np_IS.M)-(Ob_IS.F-Np_IS.F)')
paths = glue("{Dir}/Plots/GSE29718_VAT_{contrasts}.svg")
RX_grid_img(paths)
## [1] 1
## [1] 3